package ua.lviv.franko.solvers.stressanalysis.twoDimention.combined;

import java.util.ArrayList;
import java.util.Collections;
import java.util.HashSet;
import java.util.List;
import java.util.Set;

import org.jfree.data.xy.XYSeries;
import org.jfree.data.xy.XYSeriesCollection;
import org.jzy3d.maths.Range;
import org.jzy3d.plot3d.builder.Mapper;

import ua.lviv.franko.FEMShapeFunction1D.FEMShapeFunction1D1;
import ua.lviv.franko.integration.ConstantFunction1D;
import ua.lviv.franko.integration.ConstantFunction2D;
import ua.lviv.franko.integration.Gauss2D;
import ua.lviv.franko.integration.GaussQuadrature;
import ua.lviv.franko.integration.IFunction1D;
import ua.lviv.franko.integration.IFunction2D;
import ua.lviv.franko.integration.MatrixFunction1DReqtangular;
import ua.lviv.franko.integration.VectorFunction1D;
import ua.lviv.franko.matrix.IRow;
import ua.lviv.franko.matrix.RowType;
import ua.lviv.franko.matrix.hash.lu.HashMatrixLUStyle;
import ua.lviv.franko.matrix.hash.lu.LUFactorization;
import ua.lviv.franko.matrix.hash.lu.MatrixLUStyle;
import ua.lviv.franko.solvers.stressanalysis.oneDimention.AbstractElement1D;
import ua.lviv.franko.solvers.stressanalysis.oneDimention.EquationParameters;
import ua.lviv.franko.solvers.stressanalysis.oneDimention.FEElementLinear;
import ua.lviv.franko.solvers.stressanalysis.oneDimention.TimoshenkoMembraneSolver;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.Solver;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.fem.AbstractStressElement;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.fem.LinearRectangularElement;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.savestructures.bounderyconditions.BounderyConditionLinear;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.savestructures.bounderyconditions.BounderyConditionStressInterface;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.savestructures.bounderyconditions.BounderyConditionType;
import ua.lviv.franko.triangulation.SaveStructures.Node;
import ua.lviv.franko.triangulation.SaveStructures.NodePair;
import ua.lviv.franko.triangulation.SaveStructures.Rectangle;
import ua.lviv.franko.visualisation.NewOkCancelDialog;
import ua.lviv.franko.visualisation.Visualise3D;

public class MultiDimentionSolver {
	public static final int					QuadraturePointsNumber	= 7;
	protected GaussQuadrature				gaussQuadrature			= new GaussQuadrature(QuadraturePointsNumber);

	protected TimoshenkoMembraneSolver		membraneSolver;
	protected Solver						solver;
	protected List<AbstractStressElement>	nearBounderyElements;
	protected List<FEMShapeFunction1D1>		lagrangeConstraintsFunctions;

	private static final double				EPS						= 0.000000000001;
	protected HashMatrixLUStyle					mainMatrix;
	protected double[]						b;
	protected double						yValueOfTop				= 0.0;
	protected double						left					= 0.0;
	protected double						right					= 1.0;

	public MultiDimentionSolver(TimoshenkoMembraneSolver membraneSolver, Solver solver, List<AbstractStressElement> nearBounderyElements,
			List<FEMShapeFunction1D1> lagrangeConstraintsFunctions, double yValueOfTop) {
		super();
		this.membraneSolver = membraneSolver;
		this.solver = solver;
		this.nearBounderyElements = nearBounderyElements;
		this.yValueOfTop = yValueOfTop;
		this.lagrangeConstraintsFunctions = lagrangeConstraintsFunctions;
	}

	public void generate() {
		solver.generate();
		solver.applyBounderyConditions(solver.getBounderyConditions());
		membraneSolver.generate();
		membraneSolver.setBothDirihlet(new double[] { 0.0, 0.0, 0.0 }, new double[] { 0.0, 0.0, 0.0 });
		mainMatrix = new HashMatrixLUStyle(RowType.ARRAY_TYPE, solver.getB().length + membraneSolver.getB().length);
		List<Double> b = new ArrayList<Double>();
		List<IRow> rows = mainMatrix.getRows();
		List<IRow> tmpRows = solver.getMainMatrix().getRows();
		double[] tmpB = solver.getB();
		int shiftValue = tmpB.length;
		// main matrix copy
		for (int i = 0; i < tmpB.length; i++) {
			b.add(i, tmpB[i]);
			rows.set(i, tmpRows.get(i));
		}
		// membrane matrix copy
		tmpB = membraneSolver.getB();
		tmpRows = membraneSolver.getMainMatrix().getRows();
		for (int i = 0; i < tmpB.length; i++) {
			b.add(shiftValue + i, tmpB[i]);
			List<Integer> pos = tmpRows.get(i).getIndexes();
			for (int j = 0; j < pos.size(); j++) {
				pos.set(j, pos.get(j) + shiftValue);
			}
			rows.set(i + shiftValue, tmpRows.get(i));
		}
		// influence of main problem on membrane.
		// TODO: calculate sigmas ...
		for (final AbstractElement1D element1d : membraneSolver.getElementsList()) {
			List<IFunction1D> baseFunctionsList = element1d.getBaseFunctions();
			List<Integer> positions = element1d.getNodesNumbers();
			for (int i = 0; i < baseFunctionsList.size(); i++) {
				IFunction1D topFunc = baseFunctionsList.get(i);
				int rowPos = positions.get(i) * TimoshenkoMembraneSolver.NUMBER_OF_UNKNOWNS + shiftValue;
				if (positions.get(i) == 0) {
					continue;
				}
				for (int j = 0; j < nearBounderyElements.size(); j++) {
					AbstractStressElement elem = nearBounderyElements.get(j);
					double[][] res = calculateBodyForces(element1d.getLeft(), element1d.getRight(), elem.getMuCoef(), elem.getLjambdaCoef(), element1d.getEquationParameters().getK1(),
							element1d.getEquationParameters().getH(), topFunc, elem);
					List<Integer> colNumbers = elem.getNodesNumbers();
					for (int colNumPos = 0; colNumPos < colNumbers.size(); colNumPos++) {
						int colPos = colNumbers.get(colNumPos) * Solver.NUMBER_OF_UNKNOWNS;
						mainMatrix.addValue(rowPos, colPos, res[colNumPos * TimoshenkoMembraneSolver.NUMBER_OF_UNKNOWNS][0]);
						mainMatrix.addValue(rowPos + 1, colPos, res[colNumPos * TimoshenkoMembraneSolver.NUMBER_OF_UNKNOWNS + 1][0]);
						mainMatrix.addValue(rowPos + 2, colPos, res[colNumPos * TimoshenkoMembraneSolver.NUMBER_OF_UNKNOWNS + 2][0]);

						mainMatrix.addValue(rowPos, colPos + 1, res[colNumPos * TimoshenkoMembraneSolver.NUMBER_OF_UNKNOWNS][1]);
						mainMatrix.addValue(rowPos + 1, colPos + 1, res[colNumPos * TimoshenkoMembraneSolver.NUMBER_OF_UNKNOWNS + 1][1]);
						mainMatrix.addValue(rowPos + 2, colPos + 1, res[colNumPos * TimoshenkoMembraneSolver.NUMBER_OF_UNKNOWNS + 2][1]);
					}
				}
			}
		}
		// Lagrange constraints for displacements
		for (int i = 0; i < lagrangeConstraintsFunctions.size(); i++) {
			FEMShapeFunction1D1 topfunc = lagrangeConstraintsFunctions.get(i);
			mainMatrix.addRows(2);
			b.add(0.0);
			b.add(0.0);
			// influence of membrane
			for (final AbstractElement1D element1d : membraneSolver.getElementsList()) {
				double[] res = calculateLagrangeMultipliers(left, right, topfunc, element1d);
				List<Integer> positions = element1d.getNodesNumbers();

				int iterator = 0;
				for (Integer pos : positions) {
					double value = res[iterator];
					// u1 = -w
					mainMatrix.addValue(mainMatrix.getN() - 2, pos * TimoshenkoMembraneSolver.NUMBER_OF_UNKNOWNS + 1 + shiftValue, -value);
					mainMatrix.addValue(pos * TimoshenkoMembraneSolver.NUMBER_OF_UNKNOWNS + 1 + shiftValue, mainMatrix.getN() - 2, -value);
					// u2 = -v1 + h/2 * gama1
					// -v1
					mainMatrix.addValue(mainMatrix.getN() - 1, pos * TimoshenkoMembraneSolver.NUMBER_OF_UNKNOWNS + shiftValue, -value);
					mainMatrix.addValue(pos * TimoshenkoMembraneSolver.NUMBER_OF_UNKNOWNS + shiftValue, mainMatrix.getN() - 1, -value);
//					// gama
//					// get h/2
//					double hDivTwo = membraneSolver.getElementsList().get(0).getEquationParameters().getH() / 2.0;
//					mainMatrix.addValue(mainMatrix.getN() - 1, pos * TimoshenkoMembraneSolver.NUMBER_OF_UNKNOWNS + shiftValue + 2, -value * hDivTwo);
//					mainMatrix.addValue(pos * TimoshenkoMembraneSolver.NUMBER_OF_UNKNOWNS + shiftValue + 2, mainMatrix.getN() - 1, -value * hDivTwo);
					iterator++;
				}
			}
			// influence of 2d problem
			for (final AbstractStressElement tmpElement : nearBounderyElements) {
				double[] res = calculateLagrangeMultipliers(Math.max(left, topfunc.getXi() - topfunc.getH()), Math.min(right, topfunc.getXi() + topfunc.getH()), topfunc,
						tmpElement);
				List<Integer> positions = tmpElement.getNodesNumbers();

				int iterator = 0;
				for (Integer pos : positions) {
					double value = res[iterator];
					mainMatrix.addValue(mainMatrix.getN() - 2, pos * 2 + 1, value);
					mainMatrix.addValue(pos * 2 + 1, mainMatrix.getN() - 2, value);
					mainMatrix.addValue(mainMatrix.getN() - 1, pos * 2, value);
					mainMatrix.addValue(pos * 2, mainMatrix.getN() - 1, value);
					iterator++;
				}
			}
		}

		// apply membrane dirihlet conditions
		// left
		// mainMatrix.addRows(3);
		// mainMatrix.setValue(mainMatrix.getN() - 3, shiftValue, 1.0);
		// mainMatrix.setValue(mainMatrix.getN() - 2, shiftValue + 1, 1.0);
		// mainMatrix.setValue(mainMatrix.getN() - 1, shiftValue + 2, 1.0);
		//
		// mainMatrix.setValue(shiftValue, mainMatrix.getN() - 3, 1.0);
		// mainMatrix.setValue(shiftValue + 1, mainMatrix.getN() - 2, 1.0);
		// mainMatrix.setValue(shiftValue + 2, mainMatrix.getN() - 1, 1.0);
		//
		// // right
		// mainMatrix.addRows(3);
		// mainMatrix.setValue(mainMatrix.getN() - 3, tmpB.length + shiftValue - 3, 1.0);
		// mainMatrix.setValue(mainMatrix.getN() - 2, tmpB.length + shiftValue - 2, 1.0);
		// mainMatrix.setValue(mainMatrix.getN() - 1, tmpB.length + shiftValue - 1, 1.0);
		//
		// mainMatrix.setValue(tmpB.length + shiftValue - 3, mainMatrix.getN() - 3, 1.0);
		// mainMatrix.setValue(tmpB.length + shiftValue - 2, mainMatrix.getN() - 2, 1.0);
		// mainMatrix.setValue(tmpB.length + shiftValue - 1, mainMatrix.getN() - 1, 1.0);
		//
		// for (int i = 0; i < 6; i++) {
		// b.add(0.0);
		// }

		// copy to array for solving;
		this.b = new double[b.size()];
		for (int i = 0; i < this.b.length; i++) {
			this.b[i] = b.get(i);
		}

		// TODO: need to add code below to satisfy the Dirihlet conditions of 2D problem;
		// commented
		// // left
		// removeRow(shiftValue);
		// removeRow(shiftValue + 1);
		// removeRow(shiftValue + 2);
		// // right
		// removeRow(tmpB.length + shiftValue - 1);
		// removeRow(tmpB.length + shiftValue - 2);
		// removeRow(tmpB.length + shiftValue - 3);
		//
		// for (BounderyConditionStressInterface cond : solver.getBounderyConditions()) {
		// cond.applyBounderyCondition(mainMatrix, this.b);
		// for (NodePair pair : cond.getNodePairs()) {
		// removeRow(pair.getA().Number * Solver.NUMBER_OF_UNKNOWNS);
		// removeRow(pair.getB().Number * Solver.NUMBER_OF_UNKNOWNS);
		// }
		// }

		System.out.println("======================================================================");
		MatrixLUStyle.printMatrix(mainMatrix);
		System.out.println("======================================================================");
	}

	protected void removeRow(int rowPos) {
		List<Integer> positions = mainMatrix.getRows().get(rowPos).getIndexes();
		for (int i = 0; i < positions.size(); i++) {
			int pos = positions.get(i);
			if (pos == rowPos) {
				continue;
			}
			mainMatrix.setValue(pos, rowPos, 0.0);
		}
		positions.clear();
		mainMatrix.getRows().get(rowPos).getValues().clear();
		mainMatrix.setValue(rowPos, rowPos, 1.0);
	}

	protected double[] calculateLagrangeMultipliers(double left, double right, final FEMShapeFunction1D1 topfunc, final AbstractElement1D element1d) {
		return gaussQuadrature.integrate(Math.max(Math.max(left, element1d.getLeft()), topfunc.getXi() - topfunc.getH()),
				Math.min(Math.min(right, element1d.getRight()), topfunc.getXi() + topfunc.getH()), new VectorFunction1D() {

					@Override
					public int vectorSize() {
						return element1d.getBaseFunctions().size();
					}

					@Override
					public double[] calculateDerivaty(double x) {
						return null;
					}

					@Override
					public double[] calculate(double x) {
						double[] res = new double[vectorSize()];
						for (int j = 0; j < element1d.getBaseFunctions().size(); j++) {
							IFunction1D func = element1d.getBaseFunctions().get(j);
							res[j] = func.calculate(x) * topfunc.calculate(x);
						}
						return res;
					}
				});
	}

	protected double[] calculateLagrangeMultipliers(double left, double right, final IFunction1D topfunc, final AbstractStressElement element) {
		return gaussQuadrature.integrate(left, right, new VectorFunction1D() {

			@Override
			public int vectorSize() {
				return element.getBaseFunctions().size();
			}

			@Override
			public double[] calculateDerivaty(double x) {
				return null;
			}

			@Override
			public double[] calculate(double x) {
				double[] res = new double[vectorSize()];
				for (int j = 0; j < element.getBaseFunctions().size(); j++) {
					IFunction2D func = element.getBaseFunctions().get(j);
					res[j] = func.calculate(x, yValueOfTop) * topfunc.calculate(x);
				}
				return res;
			}
		});
	}

	protected double[][] calculateBodyForces(double left, double right, final double mu, final double ljambda, final IFunction1D k1, final double h, final IFunction1D topfunc,
			final AbstractStressElement element) {
		return gaussQuadrature.integrate(left, right, new MatrixFunction1DReqtangular() {

			@Override
			public int matrixRows() {
				return element.getBaseFunctions().size() * 3;
			}

			@Override
			public int matrixColumns() {
				return 2;
			}

			@Override
			public double[][] calculateDerivaty(double x) {
				return null;
			}

			@Override
			public double[][] calculate(double x) {
				double koef = 1 - (k1.calculate(x) * h) / 2.0;
				double[][] res = new double[matrixRows()][matrixColumns()];
				for (int j = 0; j < element.getBaseFunctions().size(); j++) {
					IFunction2D func = element.getBaseFunctions().get(j);
					double funcX = func.calculateDerivatyX(x, yValueOfTop);
					double funcY = func.calculateDerivatyY(x, yValueOfTop);
					double weightFunc = topfunc.calculate(x);
					res[j * 3][0] = koef * (mu + ljambda) * funcX * weightFunc;
					res[j * 3 + 1][0] = koef * mu * funcY * weightFunc;
					res[j * 3 + 2][0] = koef * (mu + ljambda) * funcX * weightFunc * h / 2.0;

					res[j * 3][1] = koef * (ljambda) * funcY * weightFunc;
					res[j * 3 + 1][1] = koef * mu * funcX * weightFunc;
					res[j * 3 + 2][1] = koef * (ljambda) * funcY * weightFunc * h / 2.0;
				}
				return res;
			}
		});
	}

	public double[] solve() {
		LUFactorization lu = new LUFactorization(this.mainMatrix);
		return lu.calculate(b);
	}

	/**
	 * @param args
	 */
	public static void main(String[] args) {
		IFunction2D f1 = new ConstantFunction2D(1.0);
		IFunction2D f3 = new ConstantFunction2D(1.0);
		IFunction1D k1 = new ConstantFunction1D(0);
		IFunction1D a1 = new ConstantFunction1D(1);
		IFunction1D sigma11Plus = new ConstantFunction1D(0);
		IFunction1D sigma13Minus = new ConstantFunction1D(0);
		IFunction1D sigma11Minus = new ConstantFunction1D(0);
		IFunction1D sigma33Minus = new ConstantFunction1D(0);
		IFunction1D sigma33Plus = new ConstantFunction1D(0);
		IFunction1D sigma13Plus = new ConstantFunction1D(0);
		double v = 0.37;
		double e = 80000;
		double v2 = e;
		double e2 = v;
		double h = 0.1;
		

		EquationParameters eqP = new EquationParameters(a1, k1, e, v, h, f1, f3, sigma13Plus, sigma13Minus, sigma11Plus, sigma11Minus, sigma33Plus, sigma33Minus);
		double a = 0.0;
		double b = 1.0;
		int n = 10;
		double hStep = (b - a) / (double) n;
		ArrayList<AbstractElement1D> elements = new ArrayList<AbstractElement1D>();
		for (int i = 0; i < n; i++) {
			elements.add(new FEElementLinear(eqP, i, a + i * hStep, a + (i + 1) * hStep));
			// elements.add(new BubbleElement(eqP, i, a + i * hStep, a + (i + 1) * hStep));
			// elements.add(new EnrichedLinearElement(eqP, i, a + i * hStep, a + (i+1) * hStep, 5, 1, new Gaussian(), 1.5));
		}

		TimoshenkoMembraneSolver s = new TimoshenkoMembraneSolver(elements);

		double ax = -0.0, bx = 1.0, ay = -0.0, by = 1.0;

		int n2 = 10, m = 10;

		double hx = (bx - ax) / (double) n2, hy = (by - ay) / (double) m;

		// double k12d = e / (1 - v * v);
		// double k22d = e * v / (1 - v * v);
		// double k32d = e / (1 - v * v) * ((1 - v) / 2.0);

		final List<AbstractStressElement> elems = new ArrayList<AbstractStressElement>();
		Set<Node> tmpNodesSet = new HashSet<Node>();

		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n2; j++) {
				Rectangle r = new Rectangle();
				r.A = new Node(ax + j * hx, ay, false, i * (n2 + 1) + j);
				r.B = new Node(ax + (j + 1) * hx, ay, false, i * (n2 + 1) + j + 1);
				r.C = new Node(ax + (j + 1) * hx, ay + hy, false, (i + 1) * (n2 + 1) + j + 1);
				r.D = new Node(ax + j * hx, ay + hy, false, (i + 1) * (n2 + 1) + j);
				tmpNodesSet.add(r.A);
				tmpNodesSet.add(r.B);
				tmpNodesSet.add(r.C);
				tmpNodesSet.add(r.D);

				// LinearRectangularElement elem = new LinearRectangularElement(r.A, r.B, r.C, r.D, k12d, k22d, k32d, f1, f3);
				LinearRectangularElement elem = new LinearRectangularElement(r.A, r.B, r.C, r.D, e2, v2, f1, f3);
				elems.add(elem);
			}
			ay += hy;
		}
		ay = -0.0;
		List<Node> nodes = new ArrayList<Node>();
		nodes.addAll(tmpNodesSet);
		Collections.sort(nodes);
		List<BounderyConditionStressInterface> bounderies = new ArrayList<BounderyConditionStressInterface>();
		List<NodePair> leftNodePairs = new ArrayList<NodePair>();
		List<NodePair> rightNodePairs = new ArrayList<NodePair>();
		List<NodePair> topNodePairs = new ArrayList<NodePair>();
		List<NodePair> bottomNodePairs = new ArrayList<NodePair>();
		for (AbstractStressElement element : elems) {
			List<Integer> positions = element.getPos();
			Node A = nodes.get(positions.get(0));
			Node B = nodes.get(positions.get(1));
			Node C = nodes.get(positions.get(2));
			Node D = nodes.get(positions.get(3));
			// left
			List<Node> tmpNodes = new ArrayList<Node>();
			if (Math.abs(A.x - ax) < EPS) {
				tmpNodes.add(A);
			}
			if (Math.abs(B.x - ax) < EPS) {
				tmpNodes.add(B);
			}
			if (Math.abs(C.x - ax) < EPS) {
				tmpNodes.add(C);
			}
			if (Math.abs(D.x - ax) < EPS) {
				tmpNodes.add(D);
			}

			if (tmpNodes.size() == 2) {
				leftNodePairs.add(new NodePair(tmpNodes.get(0), tmpNodes.get(1)));
			}

			// right
			tmpNodes.clear();
			if (Math.abs(A.x - bx) < EPS) {
				tmpNodes.add(A);
			}
			if (Math.abs(B.x - bx) < EPS) {
				tmpNodes.add(B);
			}
			if (Math.abs(C.x - bx) < EPS) {
				tmpNodes.add(C);
			}
			if (Math.abs(D.x - bx) < EPS) {
				tmpNodes.add(D);
			}

			if (tmpNodes.size() == 2) {
				rightNodePairs.add(new NodePair(tmpNodes.get(0), tmpNodes.get(1)));
			}

			// top
			tmpNodes.clear();
			if (Math.abs(A.y - by) < EPS) {
				tmpNodes.add(A);
			}
			if (Math.abs(B.y - by) < EPS) {
				tmpNodes.add(B);
			}
			if (Math.abs(C.y - by) < EPS) {
				tmpNodes.add(C);
			}
			if (Math.abs(D.y - by) < EPS) {
				tmpNodes.add(D);
			}

			if (tmpNodes.size() == 2) {
				topNodePairs.add(new NodePair(tmpNodes.get(0), tmpNodes.get(1)));
			}

			// bottom
			tmpNodes.clear();
			if (Math.abs(A.y - ay) < EPS) {
				tmpNodes.add(A);
			}
			if (Math.abs(B.y - ay) < EPS) {
				tmpNodes.add(B);
			}
			if (Math.abs(C.y - ay) < EPS) {
				tmpNodes.add(C);
			}
			if (Math.abs(D.y - ay) < EPS) {
				tmpNodes.add(D);
			}

			if (tmpNodes.size() == 2) {
				bottomNodePairs.add(new NodePair(tmpNodes.get(0), tmpNodes.get(1)));
			}
		}
		Gauss2D gauss2d = new Gauss2D(7);

		bounderies.add(new BounderyConditionLinear(BounderyConditionType.DIRIHLET_TYPE, new ConstantFunction2D(0.0), new ConstantFunction2D(0.0), leftNodePairs, gauss2d));

		bounderies.add(new BounderyConditionLinear(BounderyConditionType.DIRIHLET_TYPE, new ConstantFunction2D(0.0), new ConstantFunction2D(0.0), rightNodePairs, gauss2d));

		Solver solver = new Solver((n2 + 1) * (m + 1), elems, bounderies);
		List<AbstractStressElement> nearBounderyElements = new ArrayList<AbstractStressElement>();
		for (AbstractStressElement abstractStressElement : elems) {
			List<Integer> nodeNumbers = abstractStressElement.getNodesNumbers();
			for (Integer index : nodeNumbers) {
				if (Math.abs(nodes.get(index).y - by) < EPS) {
					nearBounderyElements.add(abstractStressElement);
					break;
				}
			}
		}
		// list for lagrange ...
		List<FEMShapeFunction1D1> femShapeFunction1D1s = new ArrayList<FEMShapeFunction1D1>();
		for (int i = 0; i <= n2; i++) {
			femShapeFunction1D1s.add(new FEMShapeFunction1D1(a + i * hx, hx));

		}
		MultiDimentionSolver multiDimentionSolver = new MultiDimentionSolver(s, solver, nearBounderyElements, femShapeFunction1D1s, by);
		multiDimentionSolver.generate();
		final double[] res = multiDimentionSolver.solve();
		for (int j = 0; j < res.length; j++) {
			System.out.println(res[j]);
		}
		System.out.println();
		Visualise3D.show(new Mapper() {
			@Override
			public double f(double x, double y) {
				double val = 0.0;
				int i = 0;
				while (val == 0.0 && i < elems.size()) {
					val += elems.get(i).calculateV(x, y, res);
					if (Math.abs(val)>EPS) {
						break;
					}
					i++;
				}
				return val;
			}
		}, new Range(ax, bx), new Range(ay, by));

		Visualise3D.show(new Mapper() {
			@Override
			public double f(double x, double y) {
				double val = 0.0;
				int i = 0;
				while (val == 0.0 && i < elems.size()) {
					val += elems.get(i).calculateU(x, y, res);
					if (Math.abs(val)>EPS) {
						break;
					}
					i++;
				}
				return val;
			}
		}, new Range(ax, bx), new Range(ay, by));

		XYSeries s1 = new XYSeries("result1", true, true);
		XYSeries s2 = new XYSeries("result2", true, true);
		XYSeries s3 = new XYSeries("result3", true, true);
		
		XYSeries u1 = new XYSeries("u1", true, true);
		XYSeries u2 = new XYSeries("u2", true, true);

		double[] tmpRes = new double[res.length - nodes.size() * 2];
		for (int i = 0; i < tmpRes.length; i++) {
			tmpRes[i] = res[nodes.size() * 2 + i];
		}

		int NNN = 20;
		AbstractElement1D el = null;
		for (int j = 0; j < elements.size(); j++) {
			el = elements.get(j);
			double hhh = hStep / NNN;
			for (int k = 0; k < NNN; k++) {
				double[] tmp = el.calculate(a + hStep * j + k * hhh, tmpRes);
				s1.add(a + hStep * j + k * hhh, tmp[0]);
				s2.add(a + hStep * j + k * hhh, tmp[1]);
				s3.add(a + hStep * j + k * hhh, tmp[2]);
				u1.add(a + hStep * j + k * hhh, tmp[1]);
				u2.add(a + hStep * j + k * hhh, tmp[0] + eqP.getH() / 2.0 * tmp[2]);
			}
		}
		double[] tmp = el.calculate(b, tmpRes);
		s1.add(b, tmp[0]);
		s2.add(b, tmp[1]);
		s3.add(b, tmp[2]);
		u1.add(b, tmp[1]);
		u2.add(b, tmp[0] + eqP.getH() / 2.0 * tmp[2]);

		XYSeriesCollection col = new XYSeriesCollection();
//		col.addSeries(s1);
//		col.addSeries(s2);
//		col.addSeries(s3);
		col.addSeries(u1);
		col.addSeries(u2);
		NewOkCancelDialog d = new NewOkCancelDialog(null, false, col);
		d.setLocationRelativeTo(null);
		d.setVisible(true);

	}
}
